#include "StdAfx.h"
#include "RNG.h"

//+------------------------------------------------------------------+
//	A levy distributed random number, with scale c and exponent alpha.
//	0 < alpha <= 2
double RNG::levy(double c, double alpha)
{
	double u, v, t, s;

	u = PI * (rand_open()-0.5);

	if(alpha == 1.0)							//	cauchy case
	{
		t = tan(u);
		return c * t;
	}

	do
	{
		v = exponential();
	}
	while(v==0);

	if(alpha == 2.0)							//	gaussian case
	{
		t = 2.0 * sin(u) * sqrt(v);
		return c * t;
	}


	double a = 1.0 / alpha,					//	general case
			 b = 1.0 - alpha;

	t = sin(alpha*u) / pow(cos(u),a);

	s = pow(cos(b*u) / v, b*a);

	return c * t * s;
}

//+------------------------------------------------------------------+
//	Levy distributed random numbers, with scale c and exponent alpha.
//	0 < alpha <= 2
__m128 RNG::levy(__m128 c, __m128 alpha)
{
	__m128 u, v, t, s;

	u = _ps_PI * (rand_open_ps()-_ps_0p5);

	v = exponential_ps();

	__m128 a = _mm_rcp_ps(alpha),
			 b = _ps_1 - alpha;

	t = sin_ps(alpha*u)    * _mm_rcp_ps(pow_ps(cos_ps(u),a));

	s = pow_ps(cos_ps(b*u) * _mm_rcp_ps(v), b*a);

	return c * t * s;
}
